Ultrarelativistic electron-hole pairing in graphene bilayer 
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We consider ground state of electron-hole graphene bilayer composed of two independently doped 
graphene layers when a condensate of spatially separated electron-hole pairs is formed. In the weak 
coupling regime the pairing affects only conduction band of electron-doped layer and valence band of 
hole-doped layer, thus the ground state is similar to ordinary BCS condensate. At strong coupling, 
an ultrarelativistic character of electron dynamics reveals and the bands which are remote from 
Fermi surfaces (valence band of electron-doped layer and conduction band of hole-doped layer) are 
also affected by the pairing. The analysis of instability of unpaired state shows that s-wave pairing 
with band-diagonal condensate structure, described by two gaps, is preferable. A relative phase of 
the gaps is fixed, however at weak coupling this fixation diminishes allowing gapped and soliton-like 
excitations. The coupled self-consistent gap equations for these two gaps are solved at zero tempera- 
ture in the constant-gap approximation and in the approximation of separable potential. It is shown 
that, if characteristic width of the pairing region is of the order of magnitude of chemical potential, 
then the value of the gap in the spectrum is not much difi'erent from the BCS estimation. However, 
if the pairing region is wider, then the gap value can be much larger and depends exponentially on 
its energy width. 

PACS numbers: 74.78.Na, 74.20.-z, 81.05.Uw 



I. INTRODUCTION 



Graphene, being a carbon layer of nionoatomic thick- 
ness with honeycomb lattice, was fabricated for the 
first time in 2004 P, 0]. Large carrier mobility makes 
graphene promising for possible nano-electronic applica- 
tions 0] • Due to symmetry properties of honeycomb 
lattice, graphene has unique electron properties: its va- 
lence and conduction bands touch at two inequivalent 
points of the first Brillouin zone, and near these points 
(in two "valleys" ) electron effective envelope wave func- 
tion obeys two-dimensional Dirac equation for massless 
particles with the Fermi velocity vp « lO^m/s w c/300 
playing the role of effective "speed of light" . This leads 
to a lot of peculiarities of graphene electron proper- 
ties different from properties of conventional quasi-two- 
dimensional electron gas 

The possibility to study massless charged fermions in 
graphene provides a bridge between condensed matter 
physics and quantum electrodynamics, allowing labora- 
tory investigation of relativistic effects which normally 
require extremal conditions Thus, it is of interest 
to study various coherent phases in graphene and super- 
conductivity in particular. As possible origins of intrinsic 
superconductivity, phonon- and plasmon-mediated mech- 
anisms in metal-coated graphene , electron correlations 
P and anisotropic electron-electron scattering near Van 
Hove singularity [9J were proposed. 

In electron systems at certain conditions, a conden- 
sation of paired electrons and holes as a result of their 
Coulomb attraction can occur, and this condensation, 
like superconductivity, is accompanied by the onset of 
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both an order parameter and the gap in single-particle 
excitation spectrum of the system. This problem was first 
considered in connection with excitonic insulator state 
however transitions between bands of paired parti- 
cles in this state lead to order parameter phase fixation 
and make a superfluidity impossible. On the contrary, 
in coupled quasi- two-dimensional quantum wells, spatial 
separation of electrons and holes strongly suppress the 
phase fixation and allow superfluid flows, which appear 
as persistent electric currents flowing in opposite direc- 
tions in different layers Experimental evidences of 
pair-condensed state formation were found in course of 
further research on electron- hole pairin g in such or simi- 
lar systems (see the reviews (Tsl [TEI Il6l|). 

The recent success in fabrication of a bilayer graphene 
system consisting of two parallel electrically indepen- 
dent grap hene layers separated by a nanometer-scale 
space [131, makes the problem of electron-hole pairing 
in graphene bilayer to be of high importance. Indepen- 
dent doping of two graphene layers would allow to achieve 
equal concentrations of electrons and holes in these lay- 
ers, whereas spatial interlayer separation D suppress in- 
terlayer tunneling and thus electron-hole recombination 
(moreover, band-structure properties of each layer are 
unaffected by the presence of the other one due to neg- 
ligible interlayer overlap of electron wave functions, in 
contrast to bilayer graphene with Bernal stacking [l8|). 
Such a system was considered first in [l^ in the weak 
coupling regime, while a possibility to achieve experi- 
mentally both weak- and strong coupling regimes was ar- 
gued. However, under increase of the coupling strength, 
a crossover to Bose-Einstein condensate (BEC) of local 
pairs, typical for attracting fermion systems (the BCS- 
BEC crossover i2Qj]), does not occur in graphene, since 
electrons and holes cannot form localized pairs due to ab- 
sence of gap in graphene spectrum [2l|, [2^ . Instead of a 
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gas of local pairs, an "ultrarelativistic" BCS state should 
be formed in graphene bilayer at strong coupling. In 
this state, electrons and holes are condensed as strongly 
overlapped Cooper pairs, as in ordinary BCS state, but 
at the same time an effectively ultrarelativistic nature of 
fermions plays an essential role. Unlike the weak cou- 
pling regime, where the pairing correlations affect only 
electrons and holes in narrow regions near their Fermi 
surfaces, in the ultrarelativistic regime particles from re- 
mote bands (valence band of electron-doped layer and 
conduction band of hole-doped one) are also involved in 
the pairing. Involvement of the remote bands into pair- 
ing process is the characteristic feature of "ultrarelativis- 
tic" BCS state in graphene bilayer in the strong coupling 
regime. In contrast to BEC of electron-hole pairs, "ul- 
trarelativistic" BCS state can possess rather high critical 
temperature and anomalous superfluid properties. 

Besides the paper [l^ , electron- hole pair condensation 
in graphene bilayer has been considered also in several 
other interesting works [H, [13, [H, [H, [13, [Ml . In [H, ^ 
the problem has been studied in the Hartree-Fock ap- 
proximation using unscreened Coulomb potential as a 
pairing interaction. Estimations [l^, [l^l of a tempera- 
ture at which the system turns into a superfluid state, ob- 
tained from numerical calculations, demonstrate that the 
superfluidity can survive up to room temperatures. How- 
ever, taking into account a screening of the interaction 
can substantially lower the transition temperature down 
to values of the order of IQ-^ K [H (see also Refs. '53, 'si'l 
concerning discussion of a role of self-consistent weaken- 
ing of the screening due to appearance of the gap). 

Note also the consideration [11] based upon exten- 
sion of the BCS theory on cone-shaped band structure, 
and the work [1^ studying condensation of electrons and 
holes with antiparallel spins emerging in graphene due 
to Zeeman splitting by in-plane magnetic field. Estima- 
tions of a highest possible critical temperature, made in 
the latter work taking into account the screening of the 
pairing Coulomb interaction, provide a result of the order 
of tens millikelvins. 

In discussing the question of quantitative prediction of 
the critical temperature of the superfluidity in electron- 
hole graphene bilayer, one should admit that, in contrast 
to the weak coupling regime, where applicability of the 
BCS theory [l^ is doubtless, at strong coupling it is nec- 
essary to take into account a number of factors going 
beyond BCS approach and capable to change essentially 
a value of the critical temperature. These factors are, 
firstly, multi-band character of the pairing and frequency 
dependence of the screened pairing interaction, and, sec- 
ondly, various correlation effects, missed in the mean-field 
approximation. Leaving a consideration of the role of dy- 
namical effects and correlations to a future work, in the 
present paper we focus our attention on multi-band char- 
acter of the pairing. We use mean-field approximation, 
and the Coulomb potential, statically screened accord- 
ing to the random phase approximation, is assumed as a 
pairing interaction. 



Some estimations of the role of multi-band character 
of BCS-like pairing of relativistic and ultrarelativistic 
fermions were made earlier both in connection with gen- 
eralization of Nambu- Jona-Lasinio model [13, [1E| and for 
superconductivity in graphene (3^ . [36| (for a simplified 
model of contact pairing potential in the latter case). 
In the present work, we perform analysis of the prob- 
lem of electron-hole pairing in graphene bilayer taking 
into account a finite range of screened pairing interac- 
tion and using physically realistic estimations for an ul- 
traviolet momentum cutoff in gap equations. By means 
of electron-hole transformation, the consideration made 
can be extended to a problem of electron-electron pairing 
in a single graphene layer (see also [13 )■ 

The paper is organized as follows: the Gor'kov equa- 
tions and gap equations are derived in Sec. [Til and solved 
approximately at T = in Sec. IIIII In Sec. lIVI the insta- 
bility of unpaired state at finite temperature is analyzed, 
and in Sec. [Vj the conclusions are made. 

II. DESCRIPTION OF THE PAIRED STATE 

Near two inequivalent points of the Brillouin zone (two 
Dirac points) the dynamics of free electrons in graphene 
is described by the equation 

VF{a ■ p}i' = Eip, (1) 

where a = {a^, ay} is a two-dimensional vector composed 
from Pauli matrices, "0 is a two-component electron wave 
function, and its components are physically the envelope 
wave functions, i.e. the amplitudes of Bloch functions 
placed on two sublattices of graphene lattice |38|. The 
solutions of ^ are the eigenfunctions ipp-yi''^ = exp(ip- 
rjfp-y having a form of a plane wave multiplied by a spinor 
/p7 — {exp(— ?y'/2),7exp(«(/j/2)}^/V2, where if is an 
azimuthal angle of the momentum vector p, measured 
from some fixed direction, 7 = ±1 is an index of a band, 
populated by the electron. The electron is located in 
the conduction band and its energy (measured from the 
Dirac point) is E — v-p\p\ at 7 +1, and in the valence 
band with E = — wf|p| at 7 = — 1. 

To describe an ultrarelativistic condensate, one could 
use a diagrammatic technique of quantum electrody- 
namics [301 in the limit of massless electrons, but a 
modification of non-relativistic diagrammatic technique 
[ig] will be used instead, since it reveals a physics of 
the multi-band pairing more clearly. For this purpose 
we introduce the following Matsubara Green functions: 
G^A^sMr) = -(T4;;,^(r)4^;+ (0)), where T is an 
ordering operator in imaginary time r, (. . .) is an av- 
erage over the thermodynamic ensemble; c^Sj^ — a^^^^, 

cSll = where ai"^^ is the destruction operator 

for electron with a momentum p, located in the layer 
a {a — 1 and 2 mean the electron- and hole-doped 
layers respectively), having a spin projection and val- 
ley denoted jointly by the index s, and populating the 
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FIG. 1: Diagrammatic representation of the Gor'kov ([Sjl and 
self-consistency Q equations. 



band 7. The switching from a to c operators corre- 
sponds to interchanging of valence and conduction bands 
in the hole-doped layer. This allows the bands of both 
kinds of pairing particles, containing their Fermi sm'- 
faces, to be denoted by 7 = +1 and named as "con- 
duction bands" , while the bands remote from the Fermi 
surfaces can be denoted by 7 = — 1 and named as "va- 
lence bands" . The bare Green functions are diagonal over 
all indices and equal to G^^^"\p, iujn) = (iw„ — 



where ^Ijy-' 



. ho 



7Df|p1 — M 3.re the energies 
of electrons and lioles measured from their chemical po- 
tentials (due to electron-hole symmetry of graphene, the 
Fermi surfaces for electrons and holes coincide at oppo- 
site chemical potentials in two layers). 

Thus, two ends of the line, denoting the Green function 



G 



(a/3) 

71S172S2 



{p,iujn), are attributed by layer (a,/?), band 



(71772) and spin-valley (si,S2) indices. Since an intra- 
and interlayer Coulomb interactions are different, we also 
assign the layer indices to each end of the interaction line 
yiaP) ^ As a result of folding over components of spinor 
electron wave functions, an additional factor, depending 
on bands and momentum directions of ingoing and out- 
going electrons, should be added to each diagram vertex. 
If the electron with the momentum p and band index 7 
enters the vertex, and the electron with the momentum 
p' and band index 7' exits, then we attribute the factor 

{fpilfpY) ^21^ 2 2 J ' 

to it, where tp and ip' are the azimuthal angles of p and 
p' respectively. Calculating an analytical expression cor- 
responding to a diagram, one should perform summation 
on layer, band and spin-valley indices in all internal ver- 
tices, taking into account that layer and spin-valley in- 
dices are conserved in the vertices, while band indices are 
not. 

Being guided by the above-stated Feynman rules, we 



draw the diagrams expressing the Gor'kov equations [4l| , 
that describe the pairing in the mean-field approxima- 
tion, and also the self-consistency conditions for the self- 
energies (Fig. [T|). The corresponding analytical expres- 
sions are 



where V'-"^\q) are the Fourier transforms of a statically 
screened intra- and interlayer electron-electron interac- 
tion. Here we do not take into account frequency de- 
pendence of self-energies, assuming that the role of dy- 
namical effects is reduced to some effective restriction of 
the pairing region in the momentum space, analogously 
to restriction of superconducting pairing in metals to the 
Debye frequency-sized region near the Fermi surface [lol | . 
Since we are interested in the interlayer pairing only, then 
only the terms of ([3]) with the anomalous self-energies 
and S(2i) are essential. The terms with the normal 
self-energies E^""^ leading to a renormalization of the 
Fermi velocity and chemical potential, will be omitted. 

The electrons and holes are spatially separated, there- 
fore particles with any spin directions and from any val- 
leys are equally suitable to be paired. This leads to 
an arbitrariness of a structure of the anomalous Green 

(12) (12) 

functions G\-^s['y2S2 and 0-^13172321 playing the role of a 
condensate wave function, on the indices si_2- In spa- 
tially homogeneous case, this spin- valley structure fac- 
torizes in a form of some (4 x 4) matrix, denoted by P: 

r-(12) _ p ^(12) ^(21) _ p+ ^(21) . 

""71S172S2 — SiS2""7l72 I ^7l*l72S2 ^ -'S1S2 Tl72- ^C 

cording to (|4]), the same factorization takes place for 



the anomalous self-energies: '^^jjsi 



^7^isi72S2 = -^sts2^7i72' Note, that the matrix P should 
be unitary — if not, then, as can be shown, excitations 
inside the gap appear and the energy of the system in- 
creases in comparison with the unitary case (thus, the 
unitarity of P means the absence of unpaired particles; 
see also (42j in connection with unitarity of a conden- 
sate). The normal Green functions are diagonal on spins 

and valleys: G^"}-y2S2 



71S172S2 



P A 

^siS2 ^7172 ^ 



"siS2'-^7l72 ■ 



Since the order parameter has a form of ?7(4) — 
U{1) X SU{4:) unitary matrix P, fluctuations of both the 
common phase of the condensate (the U{1) sector) and of 
its spin- valley structure (the SU{4) sector) are possible 
and characterized by the corresponding phase stiffnesses. 
Moreover, such topological excitations of the order pa- 
rameter, as vortices, half-vortices and quarter-vortices 
can arise (similarly to [33j . where the order parameter 
has a form of U{2) = U{1) x SU{2) matrix). 
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For simplicity, we consider the case of the gap, isotropic 
by an absolute value, i.e. | A^^^^ (p) \ does not depend on a 

direction of p. At Z-wave pairing, we have G^^^^^ (P: ^^n) ~ 
e''-'^F^^^^{p,iuj.n), G'^,^^^{p,iujn) = e~'''^F+^^ (p, iLj„) and 
A-yj^2(p) = e*''^A^j^2 (p), where (p is the azimuthal angle 

of p; the normal Green functions G^"l are independent 
on this angle. Using the aforecited notations, write down 
the systems of Gor'kov equations describing the pair- 
ing: 



(11) 

7172 



7172 
(22) 



- A^j^i^+ — ^71 72, 

A+ - 

7l7""772 



0, 



SP7lJ^' 
SP7lJ-'^^7i72 



7172 



- A+^F-y-yj — ^71 72, 
^7l7""772 



0. 



(5) 
(6) 



Performing integration on the azimuthal angle of p' in (jj]) 
and taking into account ^ , we get the self-consistent gap 
equations 

A t\ rj. ST' [ P' dp' \ I + I1I21W2 



'^,^7l72 



xVi{p,p') + ^Vi+iip,p') 

+^^-ib,p')}j^7;7^(p''*^"), (7) 



with /-wave harmonics of the pairing potential defined: 



Viip,p') = 



dip 
2^ 



V 



pi _|_ pi. 



2pp' cos ip ) e 



-ilip 



(8) 



The mixing of harmonics of different multipolarity is a 
manifestation of the "chirality" of graphene electrons @ . 

Generally, a matrix structure of the equations ([5])- 
^ is self-consistent only for certain forms of the ma- 
trix structure of the order parameter A-y-^^^. All pos- 
sible structures of the order parameter, allowing self- 
consistent solution for a contact pairing potential, were 
found in [35, 36]. It was shown, that only scalar and 
pseudoscalar order parameters has a physical meaning 
of an isotropic-gap pairing turning into usual BCS-like 
pairing in the non-relativistic limit. In our notations, 
this correspond to a gap matrix, proportional to the unit 
matrix: A. 



7172 



(X <5-yj-y2- However, the distinction of our 



finite-range pairing potential ([8]) from the contact one 
Vi{p,p') = VqSio leads to a necessity to modify this so- 
lution. We have to assume that the diagonal elements 
A^.^ can differ (similarly to [l^l). The assumption of 
the band-diagonal pairing is also the simplest physically, 
since it corresponds to a pairing of electrons and holes 
from the bands containing their Fermi surfaces with each 
other and, at the same time, to a mutual pairing of elec- 
trons and holes from the remote bands. Using V-i = Vi, 
we see from ^ and ([7]) that the band-diagonal pairing is 
self-consistent only in the case of s-wave pairing, which 



will be considered henceforth. Additional arguments for 
a band-diagonal pairing will be presented in Sec. |TVl 

The energies of quasiparticle excitations in the case of 
the band-diagonal pairing are 



Ej{p) = ^(wFP-7M)' + |A77b)|^ 

and the anomalous Green functions are 

A^^(p) 



Fjy{p, iujn) 



(9) 



(10) 



Substituting ^ into ^ and summing up over Matsub- 
ara frequencies ujn = i'2n + l)7rT, we get the following 
system of gap equations: 



A..(p) = /^{k(p,p' 



„ ^++{p') 
2E+{p') 



xtanh^-,U.(,,p')^-^tanh.^-(^') 



2T 



2T 



A-(p) = I 



p'dp' 
2tt 



VbiP,p') 



,^A++(p') 



(11) 



2E+{p') 



xtanh^+K(,,p')^-^tanh^-(^') 



2T 



2E.{j>') 



2T 



where K = (Vo + Vi)/2, H = (Vb - Vi)/2. A nonzero H 

leads to a mixing of A_|__|_ and A in pTjl and hence to 

the fixation of the relative phase of these two gaps. 

The system (1111) differs in form from the multi-band 
gap equation of [325 since the latter is obtained by a sim- 
ple extension of BCS approach on cone-like band struc- 
ture and misses the electrons chirality. In the limit of the 
contact pairing interaction, when Va = Vb, this equation 
coincides with (fTTj) . 



III. GAP AT ZERO TEMPERATURE 

A. Pairing interaction 

A solution of (fTTj) is determined by a pairing poten- 
tial V{q), which can be found taking into account two 
circumstances. Firstly, a behavior of the system is gov- 
erned by two dimcnsionless parameters [l9l | : by the ratio 
of characteristic Coulomb and quantum kinetic energies 

e2 2.19 



ehvF 



and by the dimcnsionless interlayer distance d = ppD/h, 
where pp — ii/v-p is the Fermi momentum (the sys- 
tem of coupled quantum wells is also governed by the 
same parameters, but with different [45i] : the local 
pair regime is characterized by only one dimcnsionless 
density parameter [H, Sll). Secondly, a large num- 
ber TV of fermionic flavors allows using l/A'^-expansion 
and thus validates the random-phase approximation for 
a screened potential (see [1^ and references therein). 
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In graphene bilayer, at d ^ 1, the number of flavors 
N = 2 (spin projections) x 2 (valleys) x 2 (layers) = 8, 
therefore the random-phase approximation is valid even 
at large rg [1^, |3l[. The dynamically screened inter- 
layer electron-electron interaction in this approximation 
is [13, El 



V{q,L0) 



-qD 



= (12) 



where Vq = 27re^/eg is a bare Coulomb interaction, e 
is a dielectric permittivity of a medium, surrounding 
graphene layers, and n(q,a;) is a polarization operator 
of graphene (it is the same for both layers due to the 
particle- hole symmetry). 

The cutoff energy w, specifying half-width of a pair- 
ing region around the Fermi surface, can be estimated as 
a characteristic frequency of the lower branch of plasma 
oscillations in the bilayer. This frequency, found as the 
first zero of denominator of , gives the region of fre- 
quencies where the pairing potential is attractive [loj : 



\/S,rs/d, if rsd> 1, 
w = nx<^ 8rs, if rsd< 1, < 1, (13) 
2, if rsd< 1, - 1. 

In the weak coupling regime, at small or large d, we 
have w ^ 11, and the pairing does not involve particles 
from the bands, remote from the Fermi surfaces, so we 

suppose A = 0. Therefore, at weak coupling only the 

first equation of the system (jlip remains: 



A++(p) = 



p'dp' A+ 

Va{p,P ) 



2tt 



2E+{p') 



(^')tanh^.(14) 



2T 



Due to electron chirality in graphene, this equation in- 
volves the half-sum Va — (Vq + Vi)/2 of s- and p-wave 
harmonics of the pairing potential even at s-wave pairing 
instead of only the s-wave part involved in the ordinary 
BCS gap equation [loj . 

Since w <^ fi, we can suppose the gap to be constant 
if the pairing region, in the spirit of the BCS approach, 
A++(p) = A9(w - \vFp- fj,\), where 9(a;) = at x < 0, 
Q{x) = 1 at X > 0, while the potential Va{p,p') can be 
replaced by its value Va{pF,PF) on the Fermi surface. At 
T = we obtain the asymptotical solution of in the 
BCS-like form [3: 



A = 2w exp — 



(15) 



where Aa = J^Va{pF , Pf) is the dimensionless constant of 
intraband coupling, J\f — ^/2Trhvp is the density of states 
at the Fermi level. The most favorable case for the pair- 
ing is d <^ 1, when the interlayer distance is a smallest 
characteristic distance in the system; in this case, using 
the long- wave leng th asymptotics of the static polariza- 
tion operator n(g,0) w -4Af, the potential ^ 
takes a simple form: 



UV{q,0) 



8rs 



(16) 



Integrating (fT6|) on angle according to ([8]), we get 



Aa 



27rrs — 1 — Y^lSyf — 1 arccos 



4rs 



(17) 



The expressions (fT5|) and (fTT]) determine the gap estima- 
tion in one-band pairing regime at weak coupling, consid- 
ered in fl9l|. This estimation gives ~ 10~^/i, in agreement 
with [29]. However, as will be shown below, the gap in 
a multi-band pairing regime at strong coupling can be 
much larger. 

Further, for considering the multi-band pairing, we use 
two models, allowing to reduce the integral gap equations 
(fTTj) to algebraic equations. These are a constant gap 
approximation and a separable potential approximation. 



B. Constant gap approximation 

Suppose that the gap functions Ajj{p) are constant in 
the pairing region with the energy half-width w, partly 
occupying the remote bands: A+_|_(p) = A^Q(w + n — 

vfp), a (p) = A^Q{w — ^ — vfp) (assume A+ and A_ 

to be real and positive). Furthermore, replacing the po- 
tentials Va.b by their values at the Fermi level, we reduce 
the system ([Tl]) to 



A+ = AaA+A^ 
A_ = XbA+A^ 



A+ = 



XbA-A-, 
A„A_A_, 



; dx 



2^(xt1)2 + (A±/m)2' 



(18) 



(19) 



where A^ J\fVb{pF,PF) is a dimensionless constant of 
interband coupling. Similarly to (|17p , at d <C 1 we have 



Af, = — < -27rrs 



16r2 1 
— , arccos — 

^16r2 - 1 4rs 



(20) 



According to their involvement in (|18p . the constants 
Aa and Af, are called intra- and interband coupling con- 
stants respectively. By definition, their sum Ao = Aa + Afc 
is a dimensionless s-wave harmonic of the pairing poten- 
tial at the Fermi level, whereas the difference Ai = Aa — A^ 
is a p-wave harmonic. When the range of the pairing po- 
tential is small (this occurs at large Ts), the p-wave har- 
monic tends to zero and thus Aa and At, approach each 
other. Indeed, at d < 1, r-g ^ oo ([T71), ^ give Aq ~ 1/8, 
Ai ~ l/487rrs, thus Aa,6 — > 1/16. On the contrary, when 
the range of the pairing potential is large (this is the case 
of small Ts), s- and p-wave harmonics are close by the or- 
der of magnitude, thus one cannot neglect a difference 
between Aa and A;,. Indeed, at d <C 1, rg <C 1 the expres- 
sions (Ull), ^ give Aa « (rs/7r)(-l -ln2rs), Xb « r^/Tr. 
The plots of Aq, Ai, Aa and Af, as functions of at d ^ 1, 
presented in Fig. [21 confirm these observations. 



FIG. 2: The dimensionless coupling constants at d <g 1 as 
functions of Vs'. Ao (solid line), Ai (dotted line), Aa (dashed 
line), Ai, (dash-dotted line). 



FIG. 3: The conduction-band gap A+, normalized by the 
chemical potential fi, as a function of rs at T = 0, d ^ 1 and 
various ratios of w/n, indicated near corresponding curves. 



In asymptotical case of small gaps A± <^ ii, we get 
from ^M 



A, 



w — fl 



In- 



A- 



iv — II 



In 



(21) 



Then, substituting (PTjl into (fT5|) . we obtain the asymp- 
totical expressions for the gaps: 



2/i exp 



X2 _ \2 




A_ = A+Ah 



w -j^ 



2 In ■ 



2- Aa 



^ - 1 - In 



(22) 



(23) 



These expressions differ in form from the BCS result p5)) . 
If we neglect difference between intra- and interband cou- 
pling constants and set Aa = Xb (this corresponds to a 
contact pairing potential), thenjm is reduced to the 
asymptotical formula (6) from [32l |: 



A+ = 2/icxp <^ -— H 1 

Ao i-i 



(24) 



Thus, in our model the gaps A± at T ~ depend 
on d, Tg and w. A rough estimation of the pairing re- 
gion half-width at d <C 1 is u; = max(2/x, S/Ltrg), accord- 
ing to (fT5|) . However, this estimation is only by the or- 
der of magnitude, so we consider various values of w in 



vicinity of w — Sfxrs and w = 2/i for numerical calcu- 
lations. Fig. [3] shows the conduction-band gap values 
A-|_ obtained by numerical solution of the system 
as functions of at ^ 1 and various ratios of w/fi 
(the maximum value of « 2.19/e is 2.19 in the case 
when both graphene layers are suspended in vacuum). 
We see that A^, firstly, depends on w/ii approximately 
exponentially and, secondly, quickly tends to a constant 
asymptotic ai w — const and large rs ■ This can be easily 
understood using the approximate formula (j24p and the 
plot of Ao from Fig. O on increase of rg at w = const, the 
exponent of ((24)l tends to the limiting value w/n — 17. 
Therefore, at large w the gap A_|_ depends on w expo- 
nentially, in contrast to the BCS result (fT51) . where w 
appears only as a pre-exponential factor. 

Nevertheless, according to the estimations w — 

Sfirs and w = 4/i are valid only at Ts ^ 1 and rg ~ 1 
respectively, so the appropriate values of A4. are given in 
Fig. [3] by solid lines at small rg, and by dotted lines at 
large enough Tg. Thus, in the model used the ratio w/fi 
reach 2 by the order of magnitude at maximum, and the 
maximum value of A+ is (10"'' — 10~^)/i. The gap value 
in the multi-band pairing regime in the model of constant 
gaps exceeds the BCS estimation only if w //, which is 
not the case for dynamically-screened Coulomb interac- 
tion. However, as will be shown below, going beyond the 
approximation A± = const can rise the gaps considerably 
owing to integration over wide region of momentums in 
the gap equations. 

As for the ratio of valence- to conduction-band gaps, 
the difference of these gaps A-|. and A_ is as stronger, 
as considerably the constants Aa and Xb differ, as follows 
from The difference between Aa and Xb is large 

at small Tg and small at large Tg, as noted above in the 
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as m 



2 2.19 



FIG. 4: The ratio of the valence-band gap A_ to the 
conduction-band gap A+ as a function of rs at T = 0, d ^ 1 
and various ratios of w /firs'. 16r-s (dashed line), 8rs (solid 
line), 2 (dotted line). 



discussion of Fig. [51 Fig. 3] shows the ratio A_/A+ at 
T 0, d 1 and various ratios w/ fi, as functions of rg. 
It is seen, that A_/A+ depends rather weakly on w/fi 
and, in agreement with above argumentation, increases 
and tends to unity when rg increases [i^. At rg > 1, 
the ratio A_/A_|_ is very close to unity, and this means 
that the pairings in conduction and valence bands are 
of equally importance. The multi-band character of the 
pairing can lead to a number of unusual superfluid prop- 
erties of the system. Some results concerning this prob- 
lem were obtained in [28l . [soj , where properties of vortices 
and zero-energy modes inside these vortices were studied 
in the approximation of the contact pairing potential in 
multi-band regime. 

In our case, at weak coupling, when Ah ^ Aq, the fixa- 
tion of the relative phase of A+ and A_ becomes weak, 
and in the limit Af, = the phase fixation is absent, since 
the equations in the system decouple. Therefore at 
small Tg gapped and soliton-like excitations, correspond- 
ing to oscillations of the relative phase and being solu- 
tions of the sine-Gordon equation, can arise (analogously 
to ^5lJ, see also references therein). 



C. Approximation of separable potential 

Let now try to go beyond the approximation of a con- 
stant gap and take into account dependence of the func- 
tions A±{p) on momentum. For this purpose, keep only 
a separable parts of the harmonics Vo{p,p') and Vi{p^p') 
of the statically screened potential (in the same way. 



y(-P)(p p/) ^ Vl{p,PF)VliPF,p') 

' ' Vi{pf,Pf) 

r{sep)i 



(25) 



The separable part V/'""^ {p,p') is close to the original 
potential Vi{p,p') when at least one of the momcntums p 
and p' lies in a region, which is close to the Fermi surface 
and provides the major contribution to the integrals in 
the gap equations pT|) . Thus we can suppose that the 

replacement of Vi{p,p') on its separable part Vi ^'^^\p,p') 
will not lead to essential mistakes in defining the gap 
functions A±{p). A difference between Vi ^'^^\p,p') and 
Vi (p, p') can be large if the both momentums p and p' are 
far from the Fermi surface, however, as one can see, for 
the statically screened potential in the form of (|16p the 
condition Vi ^'^^\p,p') < Vi(j),p') at any momentums ap- 
plies, so the replacement of the potential on its separable 
part can lead only to underestimation of gap values. 

Substituting ^ into we see that the functions 
A±{p) can be represented in the form 



A+{p) ^ A+Vaip) + A^Vbip), 
A^{p) ^ A + Vb{p) + A^Vaip), 



where 



Va,bip) 



vo{p) ± vi{p) 



vi{p) 



VliP,PF) 

Vi{pf,Pf) 



(26) 



(27) 



(note that A-i-(pf) — A±, so the gap in the excitation 
spectrum equals approximately to A+). As a result of 
the substitutions (|25)) - (p6| . the system (fTTjl of the inte- 
gral equations is reduced to a system of two algebraic 
equations with respect to the two gap parameters A_|_ 
and A_: 



A+ = XaAA^ 
A_ = X„CA, 



XbCA^ 
- XbAA^ 



XaCA^ + XbBA^, 
XaBA_+XtCA_, ^ 



and the coefficients of these equations 



U+ + U, 



bb^ 



PF 



bb 



u- 



Vi{p)vj{p)pdp 



2a/(wfPTm)^ 



Ai(p) 



(29) 



are functions of A± themselves, as long as A±{p) in the 
integrals are defined by ([5S)) . At large p the functions 
Va,b{p) behave as rg/2Aop, so the integrals ([29]) converge 
and we do not need any momentum cutoffs. 

In the limit A± only U^^ from the all inte- 
grals ([55)1 is logarithmically singular and diverges as 
ln(/i/A+), whereas all the other integrals converge to fi- 
nite limits. For both analytical estimations and numeri- 
cal calculations, it is convenient to isolate this singular- 
ity as follows: C/^ = ln(2^/A+) -I- [/^ and analogously 
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0.01 - 




FIG. 5: The exponent a in the approximate formula A_|. = 
2fiexp{a) for the gap parameter A+ as a function of rg. Solid 
line: the result in a framework of the separable potential ap- 
proximation according to (|32|) . Dotted line: the same, but 
with neglecting the contribution of the valence band, accord- 
ing to p3[) . Dashed line: the curve a = — 2/Ao, correspond- 
ing to BCS limit (the dash-dotted line shows the asymptotic 
a = —16 of this curve at large rg). 



A = ln(2/i/A+) + A. Denoting the hmits of A, B and C 
at A± ^ as Aq, Bq and Co respectively, we can derive 
an asymptotical solution of the system (|28p at small A±: 

A+ = 2^exp|- 1- Aa(io + So)-2A,,Co 

+ (A2 - XDiAoBo - C^o)] I {^a - (Xl - Xl)Bo] } , (30) 



A_ = A_ 



Ab + (Ag - Xl)Co 
X,^{Xl-Xl)Bo- 



(31) 



Assuming for simplicity Xa ~ Xb ~ Ao/2 (validity of 
this assumption is confirmed by Fig. [2]), we get 



A, 



2/xexp 



Ar 



+ A0 + B0 + 2Co 



(32) 



The result obtained is similar to ((24)) and differs from 
the BCS result (fTS)) by a presence of additional terms in 
the exponent; these terms are positive anywhere except 
a small region at rg ^ 1. ft is interesting to note, that, 
according to (|29)) . each of these terms can be divided 
into two parts corresponding to contributions from con- 
duction and valence bands. Excluding the valence-band 
contributions U~j from the exponent, we get an expres- 
sion for a gap, corresponding to one-band model, which 
takes into account only the conduction band: 



A^°^'^ = 2^exp<j 

Ao 



2U, 



bb 



(33) 



FIG. 6: Solid line: the gap parameter A+, normalized on 
the chemical potential fi, as a function of rs calculated in the 
separable potential approximation. Dotted line: the result 
obtained in the approximation Xa = \b = Ao/2, dashed line: 
the result for A+ in the approximation A± = const at w = 
8/irs (taken from Fig. [3} ■ 



In such a way, one part of a positive contribution to the 
exponent of ([32]) comes from the valence band, whereas 
the other part comes from the conduction band and re- 
sults from smooth dependence of A±(p) and Vi:^^^^\p,p') 
on momentums in a wide region. Fig.[5]presents the expo- 
nents a in the approximate formula A+ = 2/iexp(a) for 
the gap parameter as functions of rg in three cases; tak- 
ing into account the both valence and conduction band 
contributions, as corresponds to the formula (|32p . taking 
into account only conduction band, according to ()33p . 
and the BCS limit when the exponent is a = — 2/Ao. We 
see that, at large enough rg the additional contributions 
to the exponent from the valence and conduction bands 
are numerically large, approximately equal to each other 
and grow linearly with rg. 

The corresponding results for the gap parameter A_|_ 
are presented in Fig. [S| These results are obtained by 
numerical solving of the system (|28)) in two cases: with 
taking into account a difference between the coupling 
constants Xa.b or with neglecting this difference. The 
difference between two these curves becomes unessential 
at large Tg, in accordance with decrease of Ai relative to 
Ao (see Fig. ^ . Due to dependence of both the gap func- 
tions A±{p) and potentials Vi ^'^^\p,p') on momentums, 
the gap A+ can reach large values comparable to /j,; this 
result demonstrates a striking difference from the BCS 
resuh A+ - (10"^ - 10"^)^. 

For comparison, the curve of A-|_ in the approxima- 
tion A± = const ai w = 8/xrs, analogous to one shown 
in Fig. [31 is also shown in Fig. [Sj A good agreement 
this curve demonstrates with results obtained within the 



9 



separable potential approximation can be readily under- 
stood from the expression for the statically screened 
potential: the characteristic size of a momentum region 
providing a major contribution to the gap equations, is 
SppTs, which acts similarly to restricting the pairing re- 
gion to an energy half- width of w ~ Sfirg. 

The ratio of the gap parameters A_/A+ obtained 
by numerical calculations is close both qualitatively and 
quantitatively to the results presented in Fig. |4l 

IV. INSTABILITY OF NORMAL STATE 




Studying an instability of a normal, unpaired state 
of the system allows to determine a mean-field critical 
temperature Tc of the system transition into a super- 
fluid state. This mean-field Tc provides an upper limit 
for a Kosterlitz-Thouless transition temperature in two- 
dimensional system [ssj . Consider a vertex part of inter- 
layer electron-hole scattering in the ladder approxima- 
tion, obeying the Bethe-Salpeter equation. On decrease 
of a temperature down to the critical value T^, a pole 
appears in the vertex part as a function of a total energy 
of scattered particles. The Bethe-Salpeter equation for 
electron-hole pair at rest has a diagrammatic form shown 
in Fig. [71 and the corresponding analytical expression is 



7172 7172 



{p,p',E) = V{\p~p'\){fp^,\fp^,){U^,\fp^,) 
fir/' -. t-^ 

j£^V{\p-p"\) J2 {fp'yAf,y>,i'){fp,.\fpr',;,') 



7i 72 



'^-^ 7i 72 7i72 



{p^\p\E)J2G?^^^1{p^',E/2 



lUJ, 



xG' 



f^l{p'',-E/2- 



(34) 



where F^^^^^j^^ (p, p', E) is the vertex part describing mu- 
tual scattering of the electron from the electron-doped 
layer, located in the band 71 in the initial state and in the 
band 7^ in the final state, and the electron from the hole- 
doped layer with the initial —72 and final —72 bands; p 
and p' are the relative momentums of two electrons before 
and after scattering respectively, E is the total energy of 
electrons. At the onset of instability, the pole first arise 
at = [4(1]. The frequency sums in jMl) at = 
equal 

^ 1 - ?-^f($p7i) - »f(Cp72) 
?P7l ~^ ^P72 

where n-p{E) = [exp(i?/r) + is the Fermi distribu- 
tion. For certain band indices we have 



^77(P) = — tanh^, (35) 



1 



FIG. 7: Diagrammatic representation of the Bethe-Salpeter 
equation (p4)l . 



riP) 



sinh(^/T) 



Afi cosh(^p+/2r) cosh(^p_/2T) ' 



(36) 



We seek an instability in the s-wave scattering chan- 
nel, so F^j^2^j^^(p,p') = F^j^2^/^^(p,p',£; = 0) is inde- 
pendent on directions of the vectors p and p' , and after 
angular integration the equation ([M)) takes the form: 



1 



7l727i72 P ) 

p" dp" 

2tt 



— ; Vo{p,P ) + 



xVi{p,p') + 



E 

7^72' 



7i727"72 



Voip,p") 



717" + 7272 



Vi{p,p") 



72 7i72 



{P",P')- 



(37) 



We observe immediately that 8 components T^^^^^^^^^, 
having three band indices of one sign and the last in- 
dex of the opposite sign, identically vanish. The equa- 
tions (l37l) for remaining 8 components decouple by 
pairs onto 4 systems. Let seek their solutions in a 
form of the vertex function, constant within the pair- 
ing region: F+^3^j^^(p,p') = F+^,^j^^e(u; + fi - vpp), 



-72 7172 



{P,P') 



-727172 



Q(w — II — v-pp); the pairing 



potentials, as in the gap equations (llip . can be replaced 
by their values at the Fermi level ^54] . Then the equa- 
tions (j37|) for non-zero vertex components are 



7 



A. 



(38) 



772-72 ~ ^ 7 (39) 



where A^^^j — Xa at 71 = 72 and Ah otherwise, and the 
integrals 



(tO+7/^)/DF 



A172 ~ / S^^'^^{p)p dp 

Pf J 



(40) 
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FIG. 8: Phase diagram of the system with respect to a band 
off-diagonal pairing in the variables of the pairing region half- 
width w, normalized on the chemical potential /x, and of the 
dimensionless pairing constant A. In the region of a reentrant 
transition, an instability occurs only within some temperature 
interval. 



are introduced. 

The condition of a singularity of the vertex compo- 



nents r^j^i^i, corresponding to a band-diagonal pairing, 
is equivalent to vanishing of the determinant of the sys- 
tems (1551) and reads 



(1 - A,/++)(l - A„/__) - XII++I— = 0. (41) 

Similarly, the condition of a singularity of the compo- 
nents r^^-y^'_^', corresponding to a band off-diagonal 



pairing, follows from 



(1 - A,/+_)(l - A,/_+) - XII+-I-+ = 0. (42) 

For a qualitative analysis, assume Xa = Xh = A/2, 
then dH]) is reduced to /++ + / = 2/A and ^ — to 

-|- / |_ = 2/A. As seen from PO)) . the integral /++ 

is logarithmically divergent when T — > 0, whereas all the 

other integrals , / y, I converge to finite limits. 

This means that the instability with respect to a band- 
diagonal pairing exists at arbitrarily small A, whereas the 
band off-diagonal pairing instability occurs only when the 
coupling constant exceeds some threshold value depen- 
dent on w. 

Fig. [5] shows the phase diagram of the system with re- 
spect to the band off-diagonal pairing in the variables of 
pairing region half-width w and coupling constant A. At 
fixed w and increasing A, the stable normal state of the 
system changes into a state unstable in some tempera- 
ture interval from a nonzero minimum value to a maximal 
one. Thus, on increasing T, a reentrant phase transition 



of the system into a superfluid state occurs (by analogy 
to reentrant superconductivity). At further increase of A, 
the minimal temperature restricting a region of pairing 
instability turns into zero, and the system demonstrates 
conventional single phase transition. Such an unusual be- 
havior, involving a combination of single and reentrant 
transitions, is a consequence of unconventional interband 
Pauli-blocking term (|36p suppressing pairing in some re- 
gions of phase space at finite temperatures. In contrast, 
the intraband Pauli-blocking terms psp have a conven- 
tional hyperbolic-tangent form. 

It has been shown in Sec. IIIII that the maximal value of 
intra- and interband coupling constants Aa,b is 1/16. At 
the same time, the pairing region half-width w does not 
exceed 8^x2.19 « 17.5^i. The point w/fi = 17.5, A = 1/8 
lies far in the region of stability of the normal state in 
Fig. [51 thus an electron-hole pairing in graphene bilayer 
cannot be of band off-diagonal character. Since band- 
diagonal and band off-diagonal pairings are mutually con- 
current, then it is the former, considered in Sec. Illli that 
should be settled in the system at low temperatures. 

Solving the equation (|4ip numerically, one can find the 
critical temperature Tc concerning band-diagonal pair- 
ing. At Tc <C /i, Tc <C w we can derive an asymptotical 
formula for Tc . In this case, the integrals 1^^ (|40p are 



2/.t 



In- 



ttT 



w — II 

2fi 



InW^, 



where 7 « 0.577 is the Euler constant. Substituting these 
asymptotics into (HU, we get the result 



Tr. 



— Ah 

TT 



(43) 



coinciding with the usual BCS critical temperature to 
zero-temperature gap ratio, where A+ is given by (I22p . 
We have also solved eq. (|4T|) numerically at various rg 
and w, and obtained that the ratio of to A^{T = 0) 
is close to (US]) at rather large Tc as well. 

The Kosterlitz-Thouless superfiuid transition temper- 
ature Tkt is determined by a condition of vortex pair 
dissociation 



nps{TKT) 



(44) 



where Ps{T) is superfluid density (phase stiffness) at fi- 
nite temperature, equal to p/An at T — and vanishing 
at the mean-field transition temperature Tc [Slj . The 
factor TVkt takes into account vortex fractionalization; it 
is 1 for usual vortices and 4 for half-vortices considered 
in [33[ and having twice smaller energy then the former 
ones. In graphene bilayer, quarter-vortices are possible 
as topological excitations in f7(4) order parameter, be- 
ing matrix with respect to spin projections and valleys of 
paired particles (see Sec.|ll|. These quarter- vortices have 
four times smaller energy with respect to usual vortices, 
as can be shown by analogy with [33] , and so are the most 
favorable to be created; this leads to A'^kt = 4^ = 16 for 
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our system (which differs from the result Nkt = 4 of 
[3l|). Moreover, a participation of the remote bands in 
the pairing should lead to corrections of an expression 
for Ps{T) as compared to the one-band model, although 
at weak coupling, when one-band regime effectively takes 
place, Tkt differs insignificantly from the mean-field tem- 
perature Tc- Detailed analysis of superfluid properties of 
graphene bilayer in the multi-band regime at strong cou- 
pling and determination of corresponding Tkt will be 
addressed to future work. 



V. CONCLUSIONS 

In the present work we have considered the ultrarel- 
ativistic regime of electron-hole pairing in graphene bi- 
layer, arising at strong coupling. The characteristic fea- 
ture of the ultrarelativistic pairing is involvement of par- 
ticles from both the bands including the Fermi surfaces 
and the remote bands into the pairing correlations. Such 
a multi-band pairing has been described in the mean-field 
approximation by means of the modified diagrammatic 
technique. The Gor'kov equations obtained involve the 
gap, described by a matrix with band indices. In the 
band-diagonal pairing regime, this gap is reduced to the 
two scalar gaps: one in the conduction band and the 
other, smaller by magnitude, in the valence band. A rel- 
ative phase of these gaps is fixed, while their amplitudes 
are different as larger, as larger the first angular harmonic 
of the pairing potential with respect to the zeroth one is. 
At strong coupling, the magnitudes of the gaps are com- 
parable, so a pairing within the valence band is of equal 
importance with one in the conduction band. The fixa- 
tion of the relative phase of the gaps becomes weak at 
weak coupling, and thus gapped and soliton-like excita- 
tions, corresponding to oscillation of the relative phase, 
can arise. 

We have studied an instability of an unpaired state of 
the system appearing as a singularity of the vertex part of 
electron-hole scattering obeying the Bethe-Salpeter equa- 
tion. A search for instability in various channels, corre- 
sponding to various band structures of the vertex part, 
has led to the instability conditions for band-diagonal 
(where electron and hole are situated in the same band 
both before and after scattering) and band off-diagonal 
(when the scattering participants are situated in different 
bands) pairings. The band off-diagonal instability exists 
only when the coupling constant exceeds some thresh- 
old depending on a width of the pairing region. In some 
range of coupling strengths the instability occurs only at 
temperatures in a range from nonzero minimal to some 
maximal value, however at large enough coupling con- 
stant the lower bound of the temperature region turns 
into zero. Numerical estimations show that band off- 
diagonal instability due to screened Coulomb interac- 
tion does not realize. As for an instability in the band- 
diagonal channel, it exists at arbitrarily weak coupling 
and corresponds to the band-diagonal pairing considered 



in the paper at T = 0. Determination of Kosterlitz- 
Thouless transition temperature at strong coupling re- 
quires additional study (e.g., taking into account vortex 
fractionalization and participation of remote bands) and 
will be performed in future work. 

Estimations of a gap value at T = have been made 
in two approximations: in the approximation of constant 
gaps and in the approximation of separable potential. In 
the first case we neglect a dependence of the gap func- 
tions on momentum and restrict momentum integration 
in the gap equations by some pairing region around the 
Fermi surface. In the second case we handle a separable 
part of the statically screened interaction (which leads 
to underestimation of a gap value) and do not need any 
momentum cutoff. The numerical estimations show that 
in the first case the maximum gap value is rather small 
in agreement the BCS-type estimation [2^, if the pairing 
region half-width equals to the characteristic frequency 
of plasma oscillations in the system. Deviations from the 
BCS result owing to the multi-band pairing character can 
be significant only when the pairing region half-width is 
much larger than the chemical potential (when the pair- 
ing region half-width is not large but the pairing is multi- 
band yet, then the gap can be approximately twice larger 
than in the one-band regime, as shown in [36[). However, 
calculation results in the separable potential approxima- 
tion indicate considerable deviations from the BCS result 
due to both a dependence of the gap function and the in- 
teraction potential and integration on momentum in a 
wide region. It is shown that the separable potential ap- 
proximation yields the gap values, underestimated with 
respect to solving the integral gap equations with the 
statically screened potential; at large the gap turns out 
to be comparable to the chemical potential. This result 
indicates that more refined consideration of the problem 
is necessary in the strong coupling regime. Such a consid- 
eration should include, in particular, taking into account 
retardation and damping effects in a dynamical screen- 
ing of the electron-hole interaction and going beyond the 
mean field approximation (i.e. taking into account ver- 
tex corrections and chemical potential renormalization), 
which will be considered separately. 

Nevertheless, as is well known, the dynamically scree- 
ned potential is close to the statically screened one at 
low frequencies, tends to unscreened Coulomb potential 
at high frequencies and passes through a region of repul- 
sion concerned with plasma oscillations between these 
two limits. Since the unscreened potential acts only at 
large momentums and frequencies far from the Fermi sur- 
face, where it differs insignificantly from the statically 
screened potential, than the substitution of the overall 
potential by statically screened one is suitable in the re- 
gions both very close and very far from the Fermi surface. 
A contribution from intermediate region where the role 
of plasmons is considerable should suppress pairing, and 
quantitative study of this effect requires formulation and 
solution of Eliashberg-type equations [13, ■ 
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